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The search for elementary excitations with fractional quantum numbers is a central challenge in 
modern condensed matter physics. We explore the possibility in a realistic model for several ma- 
terials, the spin-1/2 spatially anisotropic frustrated Heisenberg antiferromagnet in two dimensions. 
By restricting the Hilbert space to that expressed by exact eigenstates of the Heisenberg chain, we 
' derive an effective Schrodinger equation valid in the weak interchain-coupling regime. The dynam- 

ical spin correlations from this approach agree quantitatively with inelastic neutron measurements 
. on the triangular antiferromagnet CS2CUCI4. The spectral features in such antiferromagnets can be 

attributed to two types of excitations: descendents of one-dimensional spinons of individual chains, 
and coherently propagating "triplon" bound states of spinon pairs. We argue that triplons are 
generic features of spatially anisotropic frustrated antiferromagnets, and arise because the bound 
spinon pair lowers its kinetic energy by propagating between chains. 
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One of the most dramatic effects of strong interactions in electronic materials is the emergence of particles with 
fractional quantum numbers, for example charge e/3 Laughlin quasiparticlcs in the fractional quantum Hall cfcct, and 
spin-charge separated excitations in one dimensional (ID) quantum wires and carbon nanotubes. Indeed, fractional- 
ization is known to be quite generic in ID conductors and magnets. 0, In this case the spin excitation carrying a 
fractional quantum number, spin 1/2, is referred to as a spinon [H, H[. In contrast, in dimensions higher than one, the 
elementary excitation from a magnetically ordered states is known as a magnon and carries spin 1 p, [^, 0] ■ Neverthe- 
1 . less, fractionalization caused by strong quantum fluctuations has been repeatedly identified theoretically as a possible 
phenomena underlying unusual experimental behavior of strongly correlated materials in two and three dimensions 
and zero magnetic field, such as high-temperature superconductors, heavy fermions, and frustrated qua ntum magnets. 
In these contexts, resonating valence bond (RVB) theories [8, 9] and slave-particle approaches [Tol. [ill [l2| have been 
developed to describe fractionalization in dimensions greater than onefr3l]. However, these approaches remain largely 
unproved. Considerable effort has been devoted to the search for such exotic behaviors for decad es 114. [lal. and only 
recently, experimental indications of fractionalized particles [TH, [l?} and disordered ground states jl8l . 19l 20 , [2~H have 
. been observed in highly frustrated antiferromagnets in two dimensions (2D). 

In this paper, we consider how spinons may appear in a 2D magnet as descendents of their ID counterparts. Our 
focus is the spin-1/2 spatially anisotropic antiferromagnetic Heisenberg model defined by the following Hamiltonian: 

<N; 

\Q . H = (JS x +l,y + J[S x ,y+l + J^Sx+l.y+l + J3SX— l,v+l) ' &x,y: (1) 

where S x , y is the spin-1/2 operator at site (x,y). Here, J denotes the intrachain coupling, and J[, J' 2 and J3 are 
interchain couplings as illustrated in Fig. [TJ We take all the coupling constants positive, reflecting antiferromagnetic 

• i-H ', interactions, focusing on the frustrated situation J[ = J' 2 + J3. The main result of this paper is a systematic method 
to calculate the inelastic magnetic structure factor S(k, u>) for the full range of energy transfers with u> varying from 
essentially zero to large scales of several times J. The result is valid provided only J' a /J is not too large, and indeed 

. . . 1 reveals characteristic features of spinon excitations. 

One strong motivation to study this model comes from experiments on the material CS2CUCI4, a spin-1/2 Heisenberg 
antiferromagnet on a spatially anisotropic triangular lattice. This corresponds to Eq.([T]) with J[ = J 2 = J' and 
J3 = 0, and the measured anisotropy is J/ J' ~ 3.[22| The spectral weight in the dynamical structure factor, S{k,u), 
measured in this compound is dominated by a broad continuum, extending up to energy above 3 J, with the usually 
strong magnon peak appearing uncharacteristically insignificant. The spectral tail for some directions in momentum 
space is well- fitted by a power-law form[l6|, [l7|]. Following this observation, numerous theories have attributed the 
behavior to fractionalized excitations of exotic two dimensional critical and/or spin liquid states. [HI, [U [HI Other 
works have compared the data to anharmonic spin wave theory. Though the latter calculations reproduce the shape 
of the observed dispersion of the (broad) peaks in CS2CUCI4, a substantial phenomenological renormalization of the 
exchange parameters must be included by hand to achieve quantitative agreement . [l^. [171 l26l. [27| However, numerical 
series expansion calculations using the un-renormalized measured J, J' values properly reproduce the experimental 
peak dispersion. [28|, [2t| 
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In this paper, we argue that the spectra in CS2CUCI4 indeed reflect the presence of spinon excitations as originally 
suggested, but that these spinons are descendents of the ID excitations of the chains formed by the strong J bonds, 
and not characteristic of any exotic 2D state. A popular argument against this notion has been that the peak energy 
has substantial dispersion in the direction transverse to the chains. We show that contrary to naive expectations, 
such dispersion does appear in a quasi-lD approach. The basic physics involved is the binding of two spinons into 
a delocalized and dispersing spin-1 pair (triplon). This is driven by kinetic energy, since only a pair of spinons may 
hop between chains. The idea is a lower dimensional analogue of Anderson's interlayer tunneling mechanism of high 
temperature superconductivity, with spinon pairs replacing Cooper pairs [30j, |31[. Triplon formation leads to specific 
signatures in the structure factor which are indeed present in the data on CsoCuCU- 

The appropriateness of the ID approach is reinforced by several works. Ref. [321 ] demonstrated that it quantitatively 
reproduces most of the complex low temperature phase diagram observed in applied magnetic fields in CS2CUCI4. It 
also showed that the frustrated J' coupling is ineffective in establishing long-range order: the characteristic energy 
scale for ordering is only of order (J') /J 3 , muchsmaller than the bare J' inter-chain exchange energy. An early 
indication of this ineffectiveness appeared in Ref. (33[, in which a "decoupled" state was suggested. More recently, the 
exact diagonalization study in Ref. [34| found that correlations between spins in neighboring chains remain extremely 
weak for J' < 0.7 J. 

This suggests that the elementary excitations (spinons) of independent spin chains are a natural basis. We therefore 
project the Hamiltonian in Eq.|T|) into the subspace of eigenstates of the ID decoupled chains (3f|[3i|- Each eigenstate 
can be characterized by the number of excited spinons, which is always even for any physical state. Remarkably, 
truncating to the first non-trivial approximation of only zero- or two-spinon states reproduces the main features of 
the spectrum of such quasi-one-dimensional frustrated antiferromagnets. Note that the two-spinon approximation is 
not a low-energy one (unlike the familiar and powerful "bosonization" technique) as it includes spinons with energies 
reaching up to 7rJ/2 3> J'. This is essential for comparison with inelastic neutron scattering data which extends over 
this full range. [l7J 

The two-spinon states of a single chain are characterized by two continuous quantum numbers, which can be 
thought of cither as the momenta k x i, k x % of the individual (unbound) spinons, or equivalently, the total momentum 
k x = k x i + k X 2 and (excitation) energy e = e s [k x i) + ts(k X 2) of the pair. We use the latter notation for convenience. 
The spinon energy is given by des Cloizeaux-Pearson dispersion, e s (k x ) = (nJ/2)\sm(k x )\ The states can also 

be characterized by their total spin and S z quantum numbers. Only the triplet (s = 1) states are relevant to the 
neutron structure factor, and one may specialize without loss of generality to the S z — +1 state, which we denote 
\k x ,e) y on chain y. For the many-chain system, the unperturbed ground state and two-spinon basis states are given 
as |G.S.)o = ®y\0) y and \k x ,e,y) = \k x ,e) y ® y '^ y \fy y >, respectively. Here, |0} y denotes the ground state of the y-th 
Heisenberg chain, of length L x . 

We choose to work with eigenstates of the total 2D momentum vector k = (k x ,k y ). Such k y eigenstates are 
superpositions: |e)k = \k x ,k y ;e) = -^=J2 v G lkyV \k x ,€,y) (here L y is the number of chains). Note that, because 

the two spinons comprising any of the original basis states always live in the same chain, there is only one intrinsic 
transverse momentum k y and not two distinct spinon momenta in the y direction. Thus there is only a one parameter 
(e) set of two-spinon states for each k x ,k y . Therefore the eigenstates in this basis take the form 

|* k ) = J tfer>fc.(e)tfk(e)|e} kj (2) 



where Dk x (e) = &(u>2,u(k x ) — e)0(e — L02,i{k x )) / ^ uj 2 u (k x ) — e 2 is the density of states of the Heisenberg chain, divided 

by L x /(2n), at momentum k x and excitation energy e 38] (0 denotes the step function). It is restricted to L<J2,i(k x ) < 
£ < ^2,u{k x ), where the boundaries of the two-spinon continuum are u>2,i(k) = e s (k x ) and u)2,u(k x ) = nJ sin[k x /2]. 
The wavefunction V'k(e) defines the spread of the eigenstate amongst this continuum. The condition that l^k) is an 
eigenstate of the Hamiltonian in the 2-spinon subspace implies the Schrodinger equation: 

«/> k (e) + / deD kx (e) J'(k)At (e)A k;o (e) ^(g) = E^ k (e), (3) 



where E is the excitation energy above the ground state, and J'(k) = 2(J{ cos k y -f- J<£ cos(k x -\- k y ^ -f- t/3 cos(k x — ^y)) 
is the Fourier transform of the interchain exchange interaction. The matrix element Ak w {e) = ~^(^\^-k yl^i e )yi 
which is crucial for this study, was obtained exactly in Ref. [3!| (see Supplementary Material). 

We solved the integral equation, Eq.([3]), numerically by carefully discretizing e to obtain a complete (in the two 
spinon subspace) set of eigenfunctions ipnk (and corresponding states l^nk)) and energies £nk, with n = 1 . . . M. Here 
we typically took the number of discretized energies M to be several thousand, as large as necessary to ensure good 
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resolution. Knowing these eigenstates, we can directly evaluate their contribution to the zero temperature dynamical 
structure factor £(k, uj): 

5(k, W )= f ^'(G.S.I^W^WIG.S^^KG.S.I^I^)! 2 ^-^). (4) 

J n 

For consistency, we approximate the ground state |G.S.) by its perturbative form to first order in J'(k), though the 
linear correction term has little effect on the results. Details are given in the Supplementary Material. 

Somewhat unexpectedly, it is possible to show analytically that the structure factor obtained in this way has 
nearly the same form as found in the well-known random phase approximation (RPA). In particular, as shown in the 
Supplementary Material, when the O(J') correction to the ground state is neglected, 

of, \ = S 1D (k x ,uj) 

l ' J [l + J'(k) X [ D (k x ,cjW + [J'Oc)x'{ D (k x ,u>)]^ {> 

Here Sxz}(k x ,oj) — x' 1 ' £) (fc :r , w)/7r = Df. x (lo^A^ (lu)\ 2 is the two-spinon structure factor of a single chain 42], and 
Xi D (k x ,uj) = J Q du'SiD (k x , a/)/(u/ — uj). This nearly coincides with the RPA expression, which is obtained by 
replacing our x with the dynamic susceptibility of a single chain, x'id ^ ^ e XiD,XiD ~> I m Xi-D- R- e Xi-D differs from 
x'\d by a small contribution from uj' < 0. However, the differences between the RPA and our two-spinon result, with 
or without the ground state correction, are very small in all situations of interest - see Supplementary Material. 
We find three types of distinctive spectral features depending on the momentum, determined by the value of J'(k): 

1. J'(k) < 0: S(k,uj) has a 5-function peak below the continuous spectrum. A typical example is shown in 
Fig. [3 (a). As discussed above, this peak arises from a triplon bound state of two spinons, l^ik)- The triplon 
dispersion ut (k) is determined from the pole of |(5J) where 

l + J'(k) X ' 1D (k x ,uj t (k))=0 (6) 

and x'lr)(k x , ^t(k)) = outside the continuum. The pole appears below u>2,i because there x'id i- s positive. 
The inter-chain dispersion of the triplon is due to the fc y -dependence of J'(k). In the weak interchain-coupling 
regime, the spectral weight Z and binding energy SE = u>2,i(k x ) — Wt(k) of the peak are small, and behave as 
Z ~ | J'(k x , k y )\ and SE ~ \J'(k x , k y )\ 2 (up to logarithmic corrections). See the Supplementary Material for 
details. 

2. J'(k) > 0: The spectral weight shifts upwards, and the peak is broadened in the continuum, see Fig. [2] (b). 
A suppression of spectral weight at the lower edge of the continuum occurs due to repulsion between the two 
spinons. When J'(k x , k y ) is sufficiently large, a (5-function peak appears above the two-spinon continuum. This 
peak corresponds to an anti-bound triplon state. However, the anti-bound peak is broadened by the four-spinon 
contribution, which leads to non-zero spectral density above the two-spinon upper-boundary, u> > ui2, u [401 ]. 

3. J'(k) = 0: For such momenta, the structure factor is identical in the two-spinon approximation to that of a 
set of decoupled chains. For the frustrated situation of principle interest, where J[ = J' 2 + Jg, this condition is 
always satisfied for k x = n (but it may also be true elsewhere). 

Now, let us compare the above features with the experimental results [H|,[l7j on CS2CUCI4. The coupling constants 
are experimentally estimated as J—0. 374(5) meV, J' = J[ = ^=0. 128(5) meV, which leads to the ratio J'/J=0.34(3) 
[22|. This compound also has some very weak additional Dzyaloshinskii-Moriya and interplane interactions not 
included in our model. These have significant effects only at very low ener gies , e.g . in inducing long-range order in 
the ground state and weak incommensurability of the ordering wave vector [32L |41| . The coupling constants of these 
interactions are experimentally estimated as about 0.05 J [22| . In this paper, we neglect them for simplicity and discuss 
the physics for e nerg ies higher than about 0.1J - note that the majority of the features in the neutron scattering 
data in Refs. fill [l7j are in this higher energy regime. In the notation of Refs. [T||[i3j the Fourier component of the 
interchain couplings reads J'(k) = 4 J' cos(fc^./2) cos(fc( / /2), where k' x and k' y are the momenta corresponding to b and 
c axes in Refs. [H,[I3, respectively: k' x = k x and k' y = k x + 2k y . 

First, we discuss the large tail of S(\s.,u>) and the interpretation of the power-law behaviors observed in CS2CUCI4 
[l7j . In the present approach, a power-law behavior at the lower edge of the continuum (102,1) is obtained only 
when J'(k) = 0. There, we expect the same behavior as occurs in decoupled Heisenberg chains, i.e. 5(k, uj) oc 
y/— ln[o> — uj 2 j]/[uj — w 2 ,;] at k x ^ ir, and S(k, uj) cx \/— Touj/lj at k x = tt near the lower edge of continuum (42j. On 
a spatially anisotropic triangular lattice, J'(k) is zero on the lines of k x = n and k y = (tt — k x )/2 in momentum space, 
which correspond to the lines of k' x — n and k' y — n. The experimental result at k' x = n is given as the G scan in Ref. 
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[r^ . The comparison of S(k, to) at k' x = tt between the present result (i.e. Sio(k, to) of the Heisenberg chain) and the 
experimental data (G scan in Fig. 5 of Ref. is shown in Fig. [31 Only a single fitting parameter - for the global 
height of intensity in this plot - has been employed. For all further comparisons (below), we will employ the same 
normalization, so the remaining comparisons are parameter-free. Although the theoretical curve and experimental 
data differ somewhat at low energies due to the neglect of long-range magnetic order and the Dzyaloshinskii-Moriya 
interaction in the theory, the agreement at higher energy is quite good. 

We next turn to the dispersion relation, which we define here, in order to ease comparison with experimental data, 
by the location of the peak w(k) in S^k, u>) at each k. A comparison of our result and the experimental data (from Fig. 
3 in Ref. Fjj) is shown in Fig. |4j It should be noted that there is no fitting parameter in this plot. The asymmetry 
of the dispersion relation of the main peak with respect to k' x = tt and 3tt observed at k' y =0 and 2tt is consistently 
reproduced by the present approach (Fig. Q] (a,b)). At k' y =3n, the dispersion relation is symmetric because J'(k) 
is zero at this momentum, which is also consistent with the experimental observation (Fig. @] (c)). Despite the ID 
starting point of the approach, it explains the experimental dependence upon transverse momentum (k' y ) as well. 
Figured] (d,e) shows «S(k, to) in the perpendicular direction to k' x at k' x = —tt/2. The sign of J'(k) changes at k' y =37T. 
This causes the following change in S'(k, to): As shown in Fig. [4] (e), a bound state is formed just below the continuum 
for k' y < 3tt. On the other hand, for k' y > 3tt, the spectral weight shifts upwards, and the peak is broadened and 
absorbed into the continuum. Put simply, the lower edge of continuum (open squares in Fig. [4| (a- d)) lies below the 
peak only in the region of J'(k) > 0, and the main peak is always observed at the lowest energy of the spectrum for 
J'(k) < 0. These features are exactly in accord with the theoretical predictions. Moreover, for J'(k) < 0, the peak is 
much sharper (in fact resolution limited) than for J'(k) > 0. This is illustrated in Figs. 2] (f,g), which compare our 
theoretical predictions to scans E,F of Ref. [l?] - note the factor of 4 larger scale in Fig. 2] (f) compared to Fig. 2] (g). 

Furthermore, the asymmetry of the experimental estimate of the upper edge of continuum with respect to k' x — tt or 
k' y = 3tt is also qualitatively understood: At the momenta with </'(k) > 0, the spectral weight shifts upwards, and the 
high-energy weight becomes larger. On the other hand, in the region of J'(k) < 0, the high-energy weight decreases, 
because part of it shifts into the bound state (Figs. [2]and|4] (e)). This feature is consistent with the behavior of the 
upper edge of continuum observed in the experiment (open circles in Fig. 21 (a-d)). Namely, the peak of the dispersion 
relation of the upper edge of continuum is observed at the momentum a little shifted toward the region of J'(k) > 
from k' x — tt or k' y = 3tt. 

Our approach allows for systematic improvements by including further multi-spinon states. As a first step, we 
included the four-spinon states in the RPA approximation. This is done numerically by expressing the matrix element 
in Eq.(2| for a finite length Heisenberg chain as a product of determinants [43l 144 , |45| . The sum rule for the total 
spectral weight and the first frequency moment is satisfied by more than 99% for the length (L x — 288) considered. We 
then calculate from this Eq.© using x'id ~* R e Xi-D an d Xxd ~^ I m Xir> an d obtain the two-dimensional 5(k, u>). We 
note that the finite-size errors for L x = 288 are insignificant compared to the instrumental resolution. The resulting 
changes are small but very encouraging - the bound state in scan E has moved down a little, making agreement with 
experimental data essentially perfect (see Fig. 21(f))- We also observe that the anti-bound states, being located in the 
region of to — k space with non-zero spectral weight for 4-spinon excitations, acquire a non-zero linewidth as expected, 
but that this is small enough that they remain visible features. 

We conclude with a general discussion of our method and its ramifications. The most significant feature is the 
emergence of a spinon bound state driven by kinetic energy. Despite the superficial similarity to the more familiar 
magnon, the physics of the bound state is quite distinct. Specifically, a magnon is a Goldstone mode which emerges 
in a long-range ordered magnet as a consequence of broken symmetry In our calculations, no such broken symmetry 
is presumed. Instead, the bound state is a true s = 1 triplet excitation, and is better characterized as a triplon than 
a magnon. The same is true for the anti-bound state. In fact, the anti-bound triplon is directly analogous to the zero 
sound mode of a neutral interacting Fermi gas(46|. 

Since in most cases, weakly coupled spin chains do eventually order at low enough temperature, it is important to 
understand the validity of our scheme in this situation. For this, it is crucial that we consider frustrated inter-chain 
couplings (Jj = J' 2 + J3). In this case, the leading divergence associated with coupling neighboring chains - the strong 
tendency to Neel order at k x — tt within each chain - is removed because J'(TT,k y ) = 0. Without this condition, 
one obtains [47l |48| strong long-range Neel order which influences spectral features on the scale of O(J'). Since this 
effect is comparable to those captured by the two-spinon approximation, the latter is unjustified without frustration. 
With frustration, any fluctuation- induced order has a much smaller characteristic energy scale [32j, |4l|, |49j, and can be 
neglected compared to the shifts of excited states captured by the present approach. Of course, the presence of any 
long-range order, however weak, does modify some excitations in a qualitative manner. The triplon, when present, 
is expected to transform smoothly into a magnon as a consequence. In regions of momentum space where no bound 
state is present below the continuum, J'(k) > 0, a magnon may weakly emerge as a consequence of long-range order. 

There are numerous important directions for extensions and applications. It would be interesting to compare 
with neutron measurements of Cs2CuBr4, which is isostructural and can be modeled similarly to CS2CUCI4 but 
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with somewhat larger J' J J w 0.5[5(|, and to search for signs of the anti-bound triplon in either material. Some 
theoretical extensions would be to include three-dimensional and Dyzaloshinskii-Moriya couplings, systematically 
treat higher-spinon states, to include thermal fluctuations at T > 0, and to take into account weak long-range order. 
A very interesting different direction is to apply analogous methods to spatially anisotropic strongly interacting 
conductors, modeled by Hubbard or t-J type Hamiltonians. Given the very small arsenal of theoretical techniques 
capable of reliably obtaining intermediate energy spectra in strongly interacting systems above one dimension, further 
investigation of such methodology seems highly worthwhile. We would like to thank J. Alicea, M.P.A. Fisher and 
R. Shindou for discussions. This work is supported by the Grant-in-aid for Scientific Research (C) No. 10354143 
from MEXT, Japan (M. K.), the Petroleum Research Fund ACS PRF 43219-AC10 (O. S.), NSF grant/DMR-0457440 
(L. B.) and the Packard Foundation (L. B.). Part of this research was completed at KITP and supported in part by 
NSF under Grant No. PHY05-51164. 

Supplementary Material 
Basis 

The two-spinon states of a single chain are characterized by two continuous quantum numbers, which can be 
thought of either as the momenta k x %, k X 2 of the individual (unbound) spinons, or equivalently, the total momentum 
k x = k x \ + k X 2 and (excitation) energy e = e s (k x i) + £s(k X 2) of the pair. The spinon excitation energy is given by the 
des Cloizeaux-Pearson dispersion relation e s (k) — (nJ/2) sin[fc] [37], which is seen to describe the lower boundary of 
the two-spinon continuum, e s (fc) = u>2,i{k). The excitation energy of the spinon pair is then expressed via the total, 
k x , and relative, q x = (k x i — k X 2)/2, momenta of the pair e(k x ,q x ) = 7rJsin[fc x /2] cos[g x ]. Observe that the upper 
(lower) boundaries of the two-spinon continuum correspond to q x = (±k x /2). We find it convenient to describe the 
two-spinon state of a chain in terms of the total momentum k x and excitation energy e of the pair. The transformation 
from q x to e explains the density of states factor -Dfc x (e) appearing in ([2]). 

To derive ([3]), we evaluate the expectation value of the Hamiltonian (JXJ) in the state ^ 

<*|H|*) k = E L x L y + J deD kx (e)e\^ k Je)\ 2 + J'(k) / dedW kx (e)D k „(e)Ai x (e)A kx (e)i;lJe)^(e), (7) 

where E = J{— In 2 + 1/4) is the ground state energy per site of decoupled chains [H, |3(|, and we made use of 
the normalization condition JdeD k:c (e)\ipk x ( e )| 2 = 1 as appropriate for the state @. "Factoring out" e-integration 
/ de.D fcx (e)V>fc x (e) in © leads to Eq.© of the main text. 

The ground state to two spinon matrix clement A kx (e) represents the key technical element of our calculation 

A fcx (e) = V2{0\Sl kx Jk x ,e) y = -iV2{0\S v _ kx Jk x ,e) y . (8) 

Its absolute value squared, M(k x , e) = |^4fc x (e)| 2 , also called the singlet-to-triplet transition rate, is obtained exactly 
by an algebraic analysis based on infinite-dimensional quantum group symmetries in Ref. [39( and subsequently 
simplified in Ref. [42|, which we follow here. The matrix element A in © is obtained as a square-root of M because 
two-spinon states with different k x and/or e are orthogonal and thus the phase can be set to zero for every set (k x , e) 
independently. Note that our choice of \k Xl e k } y as an S z = +1 cigenstate of 2 spinons in y-th. chain ensures that 
<0|5* feiJ ,|M B ,e>y = 0. 
Specifically, A k (e) = exp[-I(i)/2]/\/4^, where ^ 

sin 2 [a;t/2] 
; cosh 2 



I(t) = -la - ln|isinh(7rf/4)| + s(t) , s(t) = / dx 7^77, (9) 



Iq = 0.3677... and k, e dependence comes via the parameter t 



, r 7rf. M, u (k) -w 2 ;(fc) 

cosh ^ = V e >_, [k) ■ do) 



Integration of s(t) has been performed numerically by an adaptive quadrature algorithm. 
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Discretization in the e space 

For numerical calculations, we carefully discretize energy e for every k x value, by dividing the interval u)2 jU {k x ) ~ 
^2,i(k x ) into M discrete points. The resulting discrete M x M eigenvalue problem is then solved for every value of 
momentum k = (k x , k y ). The data points in the e-space are chosen so that the distribution of them reduces to the 
exact density of states of the Heisenberg chain [38| in the continuous limit: 

A^(e) = |l. (11) 

Next, it is convenient to define rescaled two-spinon states with Kronecker delta function overlaps in place of the Dirac 
delta-function overlapping continuum states. Specifically: 



One can check that 



M //i \M 



e e = — - e e 



J(e-e') = — — -5 ee ' = 5 ee > (13) 

2M 1 ' 2M D(e) 1 ' 2M D(e)Ae ' ' V ' 

as desired. We also rescale wavefunction as cj)(e) = \/\k x \/ (2A/)-0k(e), suppressing for brevity its dependence on the 
center-of-mass momentum k (it enters the problem only as a parameter). With these definitions, equation @ takes 
on simple form 

l*kH$>(e)| £ ) M (14) 



The Schrodinger equation, in turn, takes on matrix form 

c#e) + J'(k)M^*(e)A(e)<Ke) - E^(e). (15) 

We typically took the number of discretized energies M to be several thousand, as large as necessary to ensure good 
resolution and convergence of the results. 

Dynamical structure factor S(k, ui) 
The dynamical structure factor S(k, uj) is defined by 

S(k,«) = J ^e w *<G.S.|5° k (t)5£(0)|G.S.> ) (16) 

where no sum on a — x, y or z is implied. This is calculated within the 2-spinon subspace using the obtained 
eigenstates and energies of the effective Hamiltonian as 

S(k, U ) = lj2 l(G.S.|5Z k |vI/ k )| 2 <5^ - E k ), (17) 

E k 

where summation is over excited states of the system with momentum k and energy u> = E k . In Eq. (|17p . it is clear 
that, since the ground state is a singlet, only triplet excited states with total S z = +1 can contribute to the sum. 
Note that 

l*k> =]T^]>>(eK fe ^|^, e ,y} M (18) 



and writing S_ k = -4= J2 V > e lkyV 'S_ k we have 

5 - k l*o = Y,^T,^ e y kv{y ' v ' )s -k x J k ^y) M - ( 19 ) 



Ly 

" y,y' 



7 



This can be separated into terms with y' — y and y 1 ^ y. Projecting the resulting state into the subspace containing 
only zero or two spinons per chain, one then obtains 



|T)o + |T)i, 



with components containing zero or four total spinon excitations: 



|T>o = v^g5>(e)A fe ( e )|G.S.)o, 



(20) 



(21) 



|T)i = Ys^sjM^Yt^^^-kS^ ( 22 ) 



2M L 



y^y' 



Note that in Eq. ([22|) we have indicated explicitly the total S z of the spinon pair on chains y,y', since the two chains 
have equal and opposite S z — ±1. 

We approximate the ground state by its form to first order in J'(k), in the subspace of states with zero or two 
spinons per chain: 



|G.S.) » |G.S.) + |G.S.)i, 
with as usual, in first order perturbation theory, 

1 



IG.S. 



Eq — Ho 



n'\G.s. 



(23) 



(24) 



Here TIq = H\ j; ^ is the decoupled chain Hamiltonian, and TL 1 = ri — Tia contains the interchain exchange couplings. 
The desired matrix element then has two terms: 



(G.S.|5Z k |* k > = o(G.S.|T) + ^G.S.ITK. 
The first term is simplest. From Eq. (f2Tj) . one directly obtains 



(G.S.|T) o = ^J^^0( £ )A fcx ( e ) 



2M 



(25) 



(26) 



Now consider the second term. To evaluate this explicitly, it is useful to write 



w' = E J/ (M [~ (s£ tV s. 



-k m ,y+l + \,/-t„j+l 



QZ QZ 

D k :c ,y D -k :c ,y+l 



kx,y 



(27) 



Because |T)i contains only spin S z = ±1 spinon pairs on two chains, we can restrict the consideration of |G.S.)i to 
those components with the same structure. This means that the third term inside the square brackets in Eq. (|27[) can 
be neglected, since it creates S z — two-spinon states on the chains y, y + 1. Taking only the first two terms, we have 



e + e' 



k x .y e,e' 5=±1 

We can then evaluate the overlap. One obtains 

3/2 



1 <G.S.|T) 1 = -V2(M) J'(k)J2 



(e)A kx (e)A*_ k (e')A- kx (e>) 



(28) 



(29) 



Observe that the correction diverges at k x = 7r, where the spinons have zero energy, unless J'(k x = ir) — 0. This hints 
at the instability of the system towards magnetic ordering in the case of non-frustrated inter-chain coupling. 
Combining Eqs. (|29l26[) . the perturbation-theory- improved transition rate can be expressed as 



|(G.S.|5Z k |* k )| 2 = 



\K\ 
M 



2^ <K*)Ak* (e) < 1 - J (k) 2^ 



2M e + e' 



(30) 



This and equation (jTTJ) provides way to numerically evaluate the structure factor. Fig. [5] shows that effect of the 
perturbative correction i(G.S.|T)i on the calculated structure factor is very small. 
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Connection to random-phase approximation (RPA) 



It is possible to obtain an analytic solution to the integral equation describing the two-spinon states, Eq.([3]). To do 
so, it is convenient to use the same discretization scheme as described above in Eq. ljlip to resolve the wavefunctions 
of individual states in the continuum. The Schrodinger equation (|15p is re-written as 

e0(e) + B(E)A(e) = E<j>(e). (31) 

Note that the quantity 

B(E) = ^^J2^A*(i) (32) 

€ 

is independent of e. This allows one to completely determine the e dependence of (f)(e) as 

B(E)A(e) 

9E[e) = E _ - ■ (33) 

Inserting this form into Eq. (|32p . one obtains the eigenvalue condition 

J'M\k x \ v \A(e)\ 2 = 

2M ^ E-e [ ' 

e 

This equation has two classes of solutions. There are bound (or anti-bound) states, in which the eigenvalue E is 
well-separated from the continuum of the decoupled chains. Then the denominator in Eq. (|34p remains finite as one 
takes the discretization to zero, i.e. M — > oo. On thereby obtains the bound state condition 

J'(k) JdeD(e)^^ = l. (35) 

The second class of solution is more subtle, and occurs when E is in the range of the continuum, i.e. for finite but 
large M, it is close to one of the discretized two-spinon eigenvalues of the decoupled chains, which we denote eo- 
We assume (and confirm self-consistently) that the energy of such a state can be written as E = eo + 5/M, where 5 
remains 0(1) as M — > oo. We then rewrite Eq.(j34]) as 



J'(k)|M 



2 



e - e 5/M 



2M Wl [(eo-e) 2 - (S/M) 2 (e - e) 2 - {S/Mf 



1. (36) 



The first term contains eo — e in the numerator, and the summand is locally odd about e = eo- The contribution from 
the sum in the region where |e — eo | is 0(1 /M) is therefore negligible, because the contributions from eigenvalues e on 
either side of eo cancel. Conversely, there is a non- vanishing contribution from |e — eo| of 0(1), which in the M — > oo 
limit becomes a principle part integral. Conversely, in the second sum, there is only a 5/M factor in the numerator, 
and the integrand is locally even about e = eo. In this sum, there is an O(l) contribution from the region |e — eo | 
of 0(l/M). This must be calculated explicitly, by summing over discrete e„ = e + Aen, with integer n. Here the 
level spacing is determined from Z?(eo)Ae = \k x \/(2M), and because the sum is sharply peaked we can approximate 
A(e) s» A(eo). Because the integrand decays as |e — e |~ 2 , the contribution from |e — e | of O(l) is negligible in the 
sum, and the limits on n can be extended to ±oo. 

Carrying out this sum (using 2& 2 J^T l/( ri - 2 ~ b 2 ) = 1 — nb cot(7r6)) and taking the M — * oo limit, we find the simple 
result 

7rJ'(k) J D(e )|A(e )| 2 cot(^ 7 ) + j'(k)P / de D ^ A ^ = l, (37) 

J eo - e 



where we have defined for convenience 



2D(e )S _ SE 
\k x \ ~ Ae' 



(38) 



This form guarantees | — ^ | < 1/2, so that the shift of the energy level SE is always less than half the distance to the 
nearest eigenvalue, i.e. the levels do not cross upon increasing J' . 
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Now we turn to the determination of the structure factor. Normalization of state (fl4|) fixes ^(E 1 )! according to 



\B(E)f 



E 



\A(e)\'< 



(39) 



{E-ef_ 

This observation leads us to the structure factor, which we divide into the bound state and continuum contributions: 

S(k, w) = S bs (k, u) + 5 cont (k, u). (40) 

First consider the bound (anti-bound) state contribution - we will consider only a single bound (anti-bound) state, 
since this occurs for the triangular lattice of present interest. In any case, multiple states would simply give additive 
contributions. This gives directly a delta- function peak in the structure factor, at w = Eb s , where Eb s is the bound 
state energy: 



S hs (Ku)-\B(E)\ — 



(E-e)(E-e>) 



yr<5(w — E). 



Using Eq. 



this immediately simplifies to 



S hs (k,u) = \B(E)\ 



,2M 1 



;5{u-E). 



In this case, since E — e remains non-zero as M — > oo, the sum in Eq. (|39[) can be converted to an integral: 



2M 



deD{e) 



\A(e)\' 



{E-ef 



(41) 



(42) 



(43) 



Thus we obtain the bound state delta-function contribution 



S hs {\i,u) = |[J'(k)] 2 JdeD{e) 



\A{e)\' 



S(uj -E bs ). 



{E bs - e) 

We now turn to the continuum contribution. Following the same steps as above, we find the analog of Eq. ([42 



(44) 



(k)] : 



;S(tO-E). 



(45) 



In this case, more care must be taken in evaluating B(E) from Eq. (f39|) . because the energy denominators in the sum 
become small as M — > oo. Indeed, the sum is dominated by \E — e\ of 0(1/M), and so one may as in Eq. (f3"T]) consider 
the density of states (i.e. spacing Ae) and A(e) to be approximately constant in this region. This allows one to carry 
out the sum explicitly (using X)^°oo V( n + a ) 2 = 71-2 / sin 2 (7ra)) and obtain 



\B(E)f = 



(Ae) 2 



7T 2 CSC 2 (7T7) |^(E)| 5 



(46) 



Inserting this into Eq. (|43)) . in the large M limit the sum may be converted to an integral via ^2 E -^ Eb AI5 — > J dE, 
and thereby collapse the delta-function. One then obtains 



Scont(k,w) 



1 



1 



7T 2 csc 2 (7r 7 ) D{oj)\A{oj)\ 2 [J'(k)p 



Using esc 2 (777) = cot 2 (777) + 1 and Eq. (|3"7| . we finally arrive at 



S cont (k,w) = 



with 



X '(k x ,oj)=P 



[1 + J'(k) x '(fc., + [J'(k) X "(k x , u)} 2 
,D(J)\A(J)\ 2 ,1 



and - X "(k x ,uj) = S 1D {k x ,u) = D{w)\A{w)f 



u>' — U) IT 

The standard RPA result has the same functional form as (|4"8")) but with \' , x" replaced by Rexrpa, ImXrpa via 



Rexr P a(fc K , w) 



1 



-P 



duj 



, Im xrpa(fc a; ,a; / ) 
oj' — U) 



1 



and -Imxrpa^z, w) = sga(u)Si D (k x , \lo\ 



(47) 



(48) 



(49) 



(50) 



The close similarity of expressions and (|50[) is illustrated in FigfSl 
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Triplon dispersion in anisotropic triangular lattice 

Bound (and antibound) states outside the continuum are determined by the condition 1 + J'(k)x'(fc x , u>t(k)) = 0, 
which is just where J'(k) = 4J' cos[k x /2] cos[fc y /2] for Cs 2 CuCl 4 . 

Bound state: Since x' > f° r w < ^2,1 as follows from (|4"5)) . one needs J'(k x ,k y ) < for it to appear. Using 
the asymptotic behavior of the chain structure factor near the lower edge of the two-spinon continuum, ui2,i [42j | 
Sid ~ C^/— Eajw — cj2,z]/[^2,i(w ^ w 2,z)], we find with logarithmic accuracy 

in \ /Er^V-M^i - u>] I u) 2 ,u ~ &2,l 
X(k x ,uj^ u)2,i) = V8C 1 v arctan J — '■ '-. (51) 

Here C = exp[/o/2]/\/l67r 3 . Provided that J'(k) is negative, we readily see that the triplon binding energy behaves 
as SE = ui2,i — ^t(k) oc [J'(k)] 2 , up to very weak logarithmic corrections. The triplon appears below the continuum 
and propagates along both the x and y directions, as shown in Fig. 2a. 

To compare with the bound state data from the effective Schrodinger equation ([3]), we calculate Sid and x' 
numerically. Results for the binding energy SE = cj 2 .z — obtained in these two calculations are found to essentially 
coincide with each other as shown in Fig. [Sj We also observe that in the regions where the width of the 2-spinon 
continuum becomes comparable to SE, the scaling changes to SE ~ J'(k). The change from quadratic to linear scaling 
can be understood simply from (|51|) . and simply corresponds to the situations where the argument of the cotangent 
is large (quadratic) or 0(1) (linear). We also analyze spectral weight (residue) Z of the triplon, see (|44|) . We find 
that for small J' it scales as \/8E tx J'. A comparison between the numerical and analytical results in a wider range 
of inter-chain exchange values is shown in FigJBJ 

It should be mentioned that the multiplicative logarithmic factor in x\ eq. lfoTj) . does lead to a weak "spiral" 
instability at zero frequency - i.e. in this approach one finds a bound state with negative energy at some ordering 
momentum (measured from 77) [ill]. This instability, however, is very weak. The corresponding ordering momentum 
is extremely small, fe Xl o ~ 10~ 10 [4lj], translating into similarly tiny "instability energy" ~ Jfcx,o- Moreover, this 
classical instability is overshadowed by a stronger quantum (of collinear type) ones, of the order (J') 4 / J 3 [32|. Since 
our approximation is concerned with the features of the dynamical structure factor at energies of order J' and higher, 
we are allowed to disregard these weak instabilities. 

Anti-bound states are analyzed similarly. In this case, Wt(k) = W2,«(fe x ) + <5*(k) is above the 2-spinon continuum 
where x' < 0, see (|49[) . Since x'(w) continuously decreases in magnitude to zero but retains the same (negative) sign as 
lo in increased from lo2 U to 00, the condition for the existence of an anti-bound state is simply 1 + J (k)% (k x , uJ2,u) < 0. 
The anti-bound state therefore merges into the two-spinon continuum when 1 + J^ rit (k)x'(fca;: w 2,«) = 0. As J'(k) 
is increased above this threshold, one can show that, because S(k, w) ~ Const-y/w 2u — uj for to < t02u [42], the 
anti-bound state energy scales as 5*(k) cx (J'(k) — J^ rit (k)) 2 while its spectral weight scales as y/fit, similarly to 
the bound state situation discussed above. Unlike the bound state, the anti-bound one is not a completely sharp 
excitation. This is because it takes place in the region of 4-spinon excitations, which extend from u)2,i(k x ) up to 
w 4,ti(^r) = nJ\/2(l + I cos[fcr/2]|), 40]. Hence, x" 7^ an( i the triplon lineshape is in fact Lorentian, see (|48|) . 
However, the 4-spinon spectral weight is very small in the region between u>2. u and CJ4. U boundaries (40| . and we find 
that anti-bound states remain well defined, with the height-to-width ratio well above 1. 

Being a collective excitation above the two-particle continuum, the anti-bound triplon here is very similar to the 
familiar zero-sound mode in a Fermi-liquid, e.g. 3 He. The analogy is made much more precise by focusing on the 
region near T point in the Brillouin zone, where the 2-spinon continuum collapses onto a line: uj2, u — ^2,1 ~ k x as 
k x — > 0. In this region x'(kx,to) — \k x \/(2(v s \k x \ — cj)), where v s = TtJ/2 is spinon velocity. The dispersion is found 
immediately (see Fig. 2b): 

wt(k) = v s \k x \ (l + for J'(k) > 0. (52) 

We see that the anti-bound triplon is just an "acoustic plasmon" of the spinon gas with short-range interactions. 

As anti-bound states away from the T point require strong J' for their existence, we would like to suggest that some- 
what more two-dimensional spatially anisotropic triangular antiferromagnet Cs2CuBr4 (50j seems to be a promising 
candidate for the corresponding inelastic neutron scattering study. 
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FIG. 2: Density plot of dynamical Structure factor S(k,u>) for J' 2 — J3 = J{/2 — 0.24J at (a) k y — ty and (b) k y = 0. The 
insets show the plots at k x = tt/2. 
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co [meV] 

FIG. 3: Comparison with the experimental result for dynamical structure factor S(k, ui) at k x = tt. Solid line denotes the 
two-spinon structure factor Sid(~k,u}) of a single chain with exchange J = 0.374 meV[17J. The symbols are the experimental 
data obtained by the inelastic neutron scattering experiment on Cs2CuCk, taken from the G scan of Fig. 5 in Ref. [Ijj. The 
inset shows the log-log plot. The theoretical result is fitted to the experimental data by adjusting the height with a single 
multiplication factor. The experimental data in this and the following figures are excerpted with permission from Ref. (l7| |. 
Copyright (2003) by the American Physical Society. 
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co [meV] co [meV] 

FIG. 4: a.-d. Comparison with experimental results for dispersion relations at (a) k' y =0, (b) k' y =2n, (c) k' y =Zn and (d) 
k' x ——n/2. Density plots are the present results of dynamical structure factor S(k,uj) for J{ = J2 = 0.34J, J3 = and J=0.374 
meV. Solid and open symbols denote the main peak, and the upper and lower edges of the spectrum observed by the neutron 
scattering experiment on CS2CUCI4, respectively, taken from Ref. [l7j . Graphs (a)-(d) correspond to (1), (3), (4) and (2) of Fig. 
3 in Ref. [TtJ, respectively, e. S(k, lo) at k' x — —-k/2 near the lower edge of continuum obtained by the present approach. The 
sign of J'(k) changes at k' y — 3n. f.-g. Comparison with experimental data for the line shape of S(k, ui) at (f) k' = (~n/2, 2n) 
and (g) k' = (— tt/2,Ah). Dotted lines are present results within the 2-spinon subspace multiplied by the normalization factor 
obtained by fitting the G scan in Fig[3] Solid lines are RPA result which accounts for the 4-spinon states as obtained in a 
chain of length L x = 288, see main text for the details. The numerical data in f and g are broadened by energy resolution 
AE = 0.019 meV of the spectrometer [ItJ including the isotropic magnetic form factor of Cu 2+ ions. Symbols are experimental 
data for (f ) E scan and (g) F scan of Fig. 5 in Ref. • 
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FIG. 5: Comparison of S(k,u>) at k' = (~tt/2, 4tt) between the present approximation (|17[) an d (I30|) (blue solid line), that 
without the correction to the ground state (blue dashed line), RPA as derived in (j4]|, (|48[1 and (149 [I (red line), and the standard 
RPA (I50|) (dotted red line). 
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FIG. 6: (Left y-axis) Transition rate to the bound state (Z). Pink diamonds are obtained by the present approach, and the 
green dotted line is the analytical result in Eq. (144[) . The small deviation in the large J'/ J regime is due to the correction to the 
ground state which is not included in Eq, (l44|l . (Right y-axis) Gap between the bound state and the lower edge of continuum 
(SE). Red circles are obtained by the present approach, and the blue solid line denotes the RPA result. The data shown here 
are calculated at k = (7r/4,7r) on a spatially anisotropic triangular lattice with J[ — J^(= J') and J3 = 0. 



